Thierry Phénix & Ladislas Nalborczyk
18 juin 2019
Peut-on prédire la taille d'un individu par la taille de ses parents ?
library(tidyverse)
d1 <- read.csv("Data/parents.csv")
glimpse(d1)
Observations: 40
Variables: 4
$ gender <fct> M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, …
$ height <dbl> 62.5, 64.6, 69.1, 73.9, 67.1, 64.4, 71.1, 71.0, 67.4, 69.…
$ mother <int> 66, 58, 66, 68, 64, 62, 66, 63, 64, 65, 64, 64, 62, 69, 6…
$ father <int> 70, 69, 64, 71, 68, 66, 74, 73, 62, 69, 67, 68, 72, 66, 7…
La taille de la mère a l'air plus prédictive de la taille d'un individu, et ce d'autant plus si cet individu est une femme…
d1 %>%
gather(parent, parent.height, 3:4) %>%
ggplot(aes(x = parent.height, y = height, colour = parent, fill = parent) ) +
geom_point() +
stat_smooth(method = "lm", fullrange = TRUE) +
facet_wrap(~gender) +
theme_bw(base_size = 25)
On peut fitter plusieurs modèles avec brms, et les comparer en utilisant le WAIC.
library(brms)
m1.1 <- brm(
height ~ 1 + gender,
prior = c(
prior(normal(70, 10), class = Intercept),
prior(normal(0, 10), class = b),
prior(cauchy(0, 10), class = sigma)
),
data = d1,
iter = 4000, warmup = 1000, chains = 4, cores = 4)
m1.2 <- update(
m1.1,
newdata = d1,
formula = height ~ 1 + gender + mother + father)
m1.3 <- update(
m1.1,
newdata = d1,
formula = height ~ 1 + gender + mother + father + gender:mother)
m1.4 <- update(
m1.1,
newdata = d1,
formula = height ~ 1 + gender + mother + father + gender:father)
m1.3 %>%
plot(
pars = "^b_",
combo = c("hist", "trace"), widths = c(1, 2.5),
theme = theme_bw(base_size = 20)
)
summary(m1.1, prob = 0.89)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: height ~ 1 + gender
Data: d1 (Number of observations: 40)
Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
total post-warmup samples = 12000
Population-Level Effects:
Estimate Est.Error l-89% CI u-89% CI Eff.Sample Rhat
Intercept 63.80 0.70 62.70 64.90 11622 1.00
genderM 4.17 0.98 2.59 5.76 10153 1.00
Family Specific Parameters:
Estimate Est.Error l-89% CI u-89% CI Eff.Sample Rhat
sigma 3.11 0.37 2.57 3.74 9690 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
library(broom)
tidy(m1.1, prob = 0.89)
term estimate std.error lower upper
1 b_Intercept 63.801407 0.6954508 62.696799 64.896623
2 b_genderM 4.172265 0.9794507 2.585206 5.761208
3 sigma 3.105483 0.3720946 2.574076 3.743628
4 lp__ -109.747336 1.2538574 -112.127615 -108.395307
posterior_summary(m1.1)[1:3, ] %>%
round(digits = 3)
Estimate Est.Error Q2.5 Q97.5
b_Intercept 63.801 0.695 62.429 65.174
b_genderM 4.172 0.979 2.249 6.116
sigma 3.105 0.372 2.488 3.943
m1.1 <- add_criterion(m1.1, "waic")
m1.2 <- add_criterion(m1.2, "waic")
m1.3 <- add_criterion(m1.3, "waic")
m1.4 <- add_criterion(m1.4, "waic")
w <- loo_compare(m1.1, m1.2, m1.3, m1.4, criterion = "waic")
w[,5:8] %>% round(digits = 2)
p_waic se_p_waic waic se_waic
m1.3 5.17 1.47 187.01 10.67
m1.2 4.94 1.43 187.07 10.78
m1.4 5.33 1.57 188.06 10.90
m1.1 2.69 0.66 205.70 8.49
model_weights(m1.1, m1.2, m1.3, m1.4, weights = "waic") %>%
round(digits = 2)
m1.1 m1.2 m1.3 m1.4
0.00 0.38 0.39 0.23
my_coef_tab <-
tibble(model = c("m1.1", "m1.2", "m1.3", "m1.4")) %>%
mutate(fit = purrr::map(model, get)) %>%
mutate(tidy = purrr::map(fit, tidy)) %>%
unnest(tidy) %>%
filter(term != "lp__")
my_coef_tab %>%
complete(term = distinct(., term), model) %>%
select(model, term, estimate) %>%
mutate(estimate = round(estimate, digits = 2)) %>%
spread(key = model, value = estimate)
# A tibble: 7 x 5
term m1.1 m1.2 m1.3 m1.4
<chr> <dbl> <dbl> <dbl> <dbl>
1 b_father NA 0.18 0.18 0.16
2 b_genderM 4.17 3.54 7.3 1.3
3 b_genderM:father NA NA NA 0.03
4 b_genderM:mother NA NA -0.06 NA
5 b_Intercept 63.8 15.1 13.7 16.1
6 b_mother NA 0.570 0.6 0.570
7 sigma 3.11 2.38 2.37 2.4
# plot
my_coef_tab %>%
ggplot(aes(x = model, y = estimate, ymin = lower, ymax = upper)) +
geom_hline(yintercept = 0, color = "steelblue") +
geom_pointrange(shape = 21, color = "steelblue", fill = "steelblue") +
labs(x = NULL,
y = NULL) +
coord_flip() +
theme_classic(base_size = 30) +
theme(text = element_text(family = "Courier"),
panel.grid = element_blank(),
axis.ticks.y = element_blank(),
axis.text.y = element_text(hjust = 0),
panel.background = element_rect(),
strip.background = element_rect(color = "transparent")) +
facet_wrap(~term, ncol = 1)
pp_check(m1.3, nsamples = 1e2) + theme_bw(base_size = 20)
Les données suivantes documentent le naufrage du titanic. La colonne pclass indique la classe dans laquelle chaque passager voyageait (un proxy pour le statut socio-économique), tandis que la colonne parch indique le nombre de parents et enfants à bord.
Peut-on prédire la survie d'un passager en fonction de ces informations ?
d2 <- read.csv("Data/titanic.csv")
glimpse(d2)
Observations: 539
Variables: 5
$ survival <int> 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1…
$ pclass <fct> upper, lower, upper, lower, upper, lower, upper, upper,…
$ gender <fct> male, female, female, female, male, male, male, female,…
$ age <dbl> 22, 38, 26, 35, 35, 54, 2, 27, 4, 58, 20, 39, 14, 2, 31…
$ parch <int> 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 5, 0, 1, 0, 0, 0, 1, 5…
Cette situation revient à essayer de prédire un outcome dichotomique à l'aide de prédicteurs continus et / ou catégoriels. On peut utiliser un modèle de régression logistique (cf. cours n°6).
# centering and scaling
d2 <-
d2 %>%
mutate(
pclass = ifelse(pclass == "lower", -0.5, 0.5),
gender = ifelse(gender == "female", -0.5, 0.5),
age = scale(age) %>% as.numeric,
parch = scale(parch) %>% as.numeric
)
d2 %>%
group_by(pclass, gender) %>%
summarise(p = mean(survival) ) %>%
ggplot(aes(x = as.factor(pclass), y = p, fill = as.factor(gender) ) ) +
geom_bar(position = position_dodge(0.5), stat = "identity") +
theme_bw(base_size = 25) +
xlab("class") + ylab("p(survival)")
m0 <- brm(
survival ~ 1,
family = binomial(link = "logit"),
prior = prior(normal(0, 10), class = Intercept),
data = d2,
cores = parallel::detectCores()
)
m1 <- brm(
# using the dot is equivalent to say "all predictors" (all columns)
survival ~ .,
family = binomial(link = "logit"),
prior = c(
prior(normal(0, 10), class = Intercept),
prior(normal(0, 10), class = b)
),
data = d2,
cores = parallel::detectCores()
)
m2 <- update(m1,
formula = survival ~ 1 + pclass + gender + pclass:gender)
m3 <- update(m1,
formula = survival ~ 1 + pclass + gender + pclass:gender + age)
m3 %>%
plot(
pars = "^b_",
combo = c("hist", "trace"), widths = c(1, 2.5),
theme = theme_bw(base_size = 20)
)
summary(m3)
Family: binomial
Links: mu = logit
Formula: survival ~ pclass + gender + age + pclass:gender
Data: d2 (Number of observations: 539)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 0.32 0.18 -0.01 0.71 1551 1.00
pclass -2.97 0.40 -3.81 -2.24 1463 1.00
gender -2.61 0.36 -3.39 -1.98 1560 1.00
age -0.47 0.13 -0.74 -0.23 2861 1.00
pclass:gender 2.26 0.73 0.97 3.79 1381 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
library(broom)
tidy(m3, prob = 0.89)
term estimate std.error lower upper
1 b_Intercept 0.3243055 0.1846382 0.04603331 0.6296130
2 b_pclass -2.9662498 0.3971888 -3.61622317 -2.3670721
3 b_gender -2.6135923 0.3640584 -3.24780592 -2.0731990
4 b_age -0.4748002 0.1294763 -0.68627376 -0.2724301
5 b_pclass:gender 2.2636830 0.7257193 1.19155467 3.5041731
6 lp__ -269.8048698 1.5945055 -272.82650516 -267.8832857
posterior_summary(m3)[1:3, ] %>%
round(digits = 3)
Estimate Est.Error Q2.5 Q97.5
b_Intercept 0.324 0.185 -0.009 0.711
b_pclass -2.966 0.397 -3.806 -2.244
b_gender -2.614 0.364 -3.387 -1.978
m0 <- add_criterion(m0, "waic")
m1 <- add_criterion(m1, "waic")
m2 <- add_criterion(m2, "waic")
m3 <- add_criterion(m3, "waic")
w <- loo_compare(m0, m1, m2, m3, criterion = "waic")
w[,5:8] %>% round(digits = 2)
p_waic se_p_waic waic se_waic
m3 5.33 0.80 512.69 26.52
m1 5.05 0.45 520.99 25.91
m2 4.24 0.66 524.75 25.42
m0 0.97 0.02 717.98 10.99
model_weights(m0, m1, m2, m3, weights = "waic") %>%
round(digits = 2)
m0 m1 m2 m3
0.00 0.02 0.00 0.98
pp_check(m3, nsamples = 1e2) + theme_bw(base_size = 25)
Ce jeu de données recense des informations sur le diamètre (colonne diam) de 80 pommes (chaque pomme étant identifiée par la colonne id), poussant sur 10 arbres différents (colonne tree). On a mesuré ce diamètre pendant 6 semaines successives (colonne time).
Que peut-on dire de la pousse de ces pommes, tout en considérant les structures de dépendance existant dans les données (chaque pomme poussait sur un arbre différent) ?
d3 <- read.csv("Data/apples.csv")
glimpse(d3)
Observations: 480
Variables: 5
$ tree <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ apple <int> 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 10, …
$ id <int> 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 10, …
$ time <int> 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2…
$ diam <dbl> 2.90, 2.90, 2.90, 2.93, 2.94, 2.94, 2.86, 2.90, 2.93, 2.96…
Cette situation revient à essayer de prédire un outcome dichotomique à l'aide de prédicteurs continus et / ou catégoriels. On peut utiliser un modèle multi-niveaux (ou modèle mixte, cf. cours n°8).
d3 <- d3 %>% filter(diam != 0) # removing missing data
d3 %>% # plotting
ggplot(aes(x = time, y = diam, colour = as.factor(apple) ) ) +
geom_point(show.legend = FALSE) +
geom_line(show.legend = FALSE) +
facet_wrap(~tree, ncol = 5) +
theme_bw(base_size = 25)
On peut fitter plusieurs modèles avec brms et les comparer ensuite en utilisant le WAIC.
p1 <- c(
prior(normal(0, 10), class = Intercept),
prior(cauchy(0, 10), class = sigma)
)
m1 <- brm(
diam ~ 1,
prior = p1,
data = d3,
cores = parallel::detectCores()
)
p2 <- c(
prior(normal(0, 10), class = Intercept),
prior(normal(0, 10), class = b),
prior(cauchy(0, 10), class = sigma)
)
m2 <- brm(
diam ~ 1 + time,
prior = p2,
data = d3,
cores = parallel::detectCores()
)
p3 <- c(
prior(normal(0, 10), class = Intercept),
prior(normal(0, 10), class = b),
prior(cauchy(0, 10), class = sd),
prior(cauchy(0, 10), class = sigma)
)
m3 <- brm(
diam ~ 1 + time + (1 | tree),
prior = p3,
data = d3,
cores = parallel::detectCores()
)
p4 <- c(
prior(normal(0, 10), class = Intercept),
prior(normal(0, 10), class = b),
prior(cauchy(0, 10), class = sd),
prior(cauchy(0, 10), class = sigma),
prior(lkj(2), class = cor)
)
m4 <- brm(
diam ~ 1 + time + (1 + time | tree),
prior = p4,
data = d3,
cores = parallel::detectCores(),
control = list(adapt_delta = 0.99)
)
p5 <- c(
prior(normal(0, 10), class = Intercept),
prior(normal(0, 10), class = b),
prior(cauchy(0, 10), class = sd),
prior(cauchy(0, 10), class = sigma),
prior(lkj(2), class = cor)
)
m5 <- brm(
diam ~ 1 + time + (1 + time | tree / apple),
prior = p5,
data = d3,
cores = parallel::detectCores(),
control = list(adapt_delta = 0.99)
)
m1 <- add_criterion(m1, "waic")
m2 <- add_criterion(m2, "waic")
m3 <- add_criterion(m3, "waic")
m4 <- add_criterion(m4, "waic")
m5 <- add_criterion(m5, "waic")
w <- loo_compare(m1, m2, m3, m4, m5, criterion = "waic")
w[,5:8] %>% round(digits = 2)
p_waic se_p_waic waic se_waic
m5 114.76 7.92 -2297.34 33.76
m3 12.25 1.54 -824.42 43.24
m4 13.71 1.66 -822.92 43.26
m2 4.38 1.02 -778.63 48.58
m1 2.94 0.77 -701.90 43.54
model_weights(m1, m2, m3, m4, m5, weights = "waic") %>%
round(digits = 2)
m1 m2 m3 m4 m5
0 0 0 0 1
summary(m5)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: diam ~ 1 + time + (1 + time | tree/apple)
Data: d3 (Number of observations: 451)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~tree (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 0.02 0.01 0.00 0.05 339 1.00
sd(time) 0.00 0.00 0.00 0.01 1703 1.00
cor(Intercept,time) 0.20 0.42 -0.68 0.88 708 1.00
~tree:apple (Number of levels: 80)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 0.09 0.01 0.08 0.11 846 1.00
sd(time) 0.01 0.00 0.00 0.01 1630 1.00
cor(Intercept,time) 0.50 0.13 0.23 0.73 2201 1.00
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 2.83 0.01 2.81 2.86 661 1.01
time 0.03 0.00 0.02 0.03 1902 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 0.02 0.00 0.01 0.02 2624 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
post <- posterior_samples(m5, "b")
ggplot(data = d3, aes(x = time, y = diam) ) +
geom_point(alpha = 0.5, shape = 1) +
geom_abline(
data = post, aes(intercept = b_Intercept, slope = b_time),
alpha = 0.01, size = 0.5) +
theme_bw(base_size = 25)
library(tidybayes)
library(modelr)
d3 %>%
group_by(tree, apple) %>%
data_grid(time = seq_range(time, n = 1e2) ) %>%
add_fitted_samples(m5, n = 1e2) %>%
ggplot(aes(x = time, y = diam, colour = factor(apple) ) ) +
geom_line(
aes(y = estimate, group = paste(apple, .iteration) ),
alpha = 0.2, show.legend = FALSE) +
facet_wrap(~tree, ncol = 5) +
theme_bw(base_size = 25)
pp_check(m3, nsamples = 1e2) + theme_bw(base_size = 25)
pp_check(m5, nsamples = 1e2) + theme_bw(base_size = 25)
Ces données recensent le nombre de candidatures pour 6 départements (colonne dept) à Berkeley (données disponibles dans le package rethinking). La colonne admit indique le nombre de candidatures acceptées et la colonne reject le nombre de candidatures rejetées (la colonne applications est simplement la somme des deux), en fonction du genre des candidats (applicant.gender).
On veut savoir s'il existe un biais de genre dans l'admission des étudiants à Berkeley.
d4 <- read.csv("Data/UCBadmit.csv")
glimpse(d4)
Observations: 12
Variables: 6
$ X <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12
$ dept <fct> A, A, B, B, C, C, D, D, E, E, F, F
$ applicant.gender <fct> male, female, male, female, male, female, male,…
$ admit <int> 512, 89, 353, 17, 120, 202, 138, 131, 53, 94, 2…
$ reject <int> 313, 19, 207, 8, 205, 391, 279, 244, 138, 299, …
$ applications <int> 825, 108, 560, 25, 325, 593, 417, 375, 191, 393…
Cette situation revient à essayer de prédire un outcome dichotomique (admit, reject) à l'aide de prédicteurs continus et / ou catégoriels.
d4 %>%
ggplot(aes(x = dept, y = admit / applications) ) +
geom_bar(stat = "identity") +
facet_wrap(~applicant.gender) +
theme_bw(base_size = 25)
On peut fitter plusieurs modèles avec brms et les comparer ensuite en utilisant le WAIC.
# centering gender predictor
d4$gender <- ifelse(d4$applicant.gender == "female", -0.5, 0.5)
# creating an index for department
d4$dept_id <- coerce_index(d4$dept)
p1 <- c(
prior(normal(0, 10), class = "Intercept"),
prior(cauchy(0, 2), class = "sd")
)
m1 <- brm(
admit | trials(applications) ~ 1 + (1 | dept_id),
data = d4, family = binomial,
prior = p1,
warmup = 1000, iter = 5000,
control = list(adapt_delta = 0.99, max_treedepth = 12)
)
p2 <- c(
prior(normal(0, 10), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(cauchy(0, 2), class = "sd")
)
m2 <- brm(
admit | trials(applications) ~ 1 + gender + (1 | dept_id),
data = d4, family = binomial,
prior = p2,
warmup = 1000, iter = 5000,
control = list(adapt_delta = 0.99, max_treedepth = 12)
)
p3 <- c(
prior(normal(0, 10), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(cauchy(0, 2), class = "sd"),
prior(lkj(2), class = "cor")
)
m3 <- brm(
admit | trials(applications) ~ 1 + gender + (1 + gender | dept_id),
data = d4, family = binomial,
prior = p3,
warmup = 1000, iter = 5000,
control = list(adapt_delta = 0.99, max_treedepth = 12)
)
m1 <- add_criterion(m1, "waic")
m2 <- add_criterion(m2, "waic")
m3 <- add_criterion(m3, "waic")
w <- loo_compare(m1, m2, m3, criterion = "waic")
w[,5:8] %>% round(digits = 2)
p_waic se_p_waic waic se_waic
m3 6.85 1.32 91.15 4.60
m1 6.51 2.28 105.13 17.95
m2 9.37 3.04 108.39 16.54
model_weights(m1, m2, m3, weights = "waic") %>%
round(digits = 2)
m1 m2 m3
0 0 1
stanplot(m3, pars = c("^b_", "^sd_")) +
theme_fivethirtyeight(base_size = 25) +
theme(axis.text.y = element_text(hjust = 0))
library(broom)
tidy(m3, par_type = "non-varying", prob = 0.95)
term estimate std.error lower upper
1 Intercept -0.5800742 0.6772113 -1.9556560 0.7947479
2 gender -0.1622283 0.2471726 -0.6639525 0.3396001
library(tidybayes)
library(modelr)
d4 %>%
group_by(dept_id, applications) %>%
data_grid(gender = seq_range(gender, n = 1e2) ) %>%
add_fitted_samples(m3, newdata = ., n = 100, scale = "linear") %>%
mutate(estimate = plogis(estimate) ) %>%
ggplot(aes(x = gender, y = estimate, group = .iteration) ) +
geom_hline(yintercept = 0.5, lty = 2) +
geom_line(aes(y = estimate, group = .iteration), size = 0.5, alpha = 0.1) +
facet_wrap(~dept_id, nrow = 2) +
theme_bw(base_size = 20)
pp_check(m3, nsamples = 1e2) + theme_bw(base_size = 25)
Le dilemme du tramway (trolley problem) est une expérience de pensée qui permet d'étudier les déterminants des jugements de moralité (i.e., qu'est-ce qui fait qu'on juge une action comme morale, ou pas ?).
Sous une forme générale, ce dilemme consiste à poser la question suivante: si une personne peut effectuer un geste qui bénéficiera à un groupe de personnes A, mais, ce faisant, nuira à une personne B (seule); est-il moral pour la personne d'effectuer ce geste ?
Voir ce lien pour plus d'informations.
Généralement, on fait lire des scénarios aux participants de l'étude, dans lesquels un individu doit prendre une décision dans une situation similaire à celle décrite à la slide précédente. Par exemple, imaginons que Denis ait le choix entre ne rien faire et laisser un train tuer cinq personnes, ou faire dérailler ce train mais tuer une personne. Ensuite, on demande aux participants de juger de la moralité de l'action choisie de Denis, sur une échelle de 1 à 7.
Des études antérieures ont montré que ces jugements de moralité sont grandement influencés par trois mécanismes de raisonnement inconscients:
Ce jeu de données comprend 12 colonnes et 9930 lignes, pour 331 individus. L'outcome qui nous intéresse est response, un entier pouvant aller de 1 à 7, qui indique à quel point il est permis (moralement) de réaliser l'action décrite dans le scénario correspondant, en fonction de l'age et genre (male) du participant (id).
On se demande comment les jugements d'acceptabilité sont influencés par les trois principes décrits slide précédente. Ces trois principes correspondent aux trois variables, action, intention, et contact (dummy-coded).
d5 <- read.csv("Data/morale.csv")
glimpse(d5)
Observations: 9,930
Variables: 7
$ response <int> 4, 3, 4, 3, 3, 3, 5, 4, 4, 4, 4, 4, 4, 5, 4, 4, 4, 4, …
$ id <fct> 96;434, 96;434, 96;434, 96;434, 96;434, 96;434, 96;434…
$ age <int> 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14…
$ male <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ action <int> 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, …
$ intention <int> 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ contact <int> 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
On essaye de prédire un jugement exprimé sous forme d'entier entre 1 et 7. Autrement dit, la variable qu'on essaye de prédire est une variable catégorielle, et dont les catégories sont ordonnées de 1 à 7…
d5$response %>% table %>% plot(xlab = "response", ylab = "", cex.axis = 2, cex.lab = 2)
Ce type de données peut se modéliser en utilisant le modèle de régression logistique ordinale. Ci-dessous un exemple en utilisant brms, et les priors par défaut (NB: ces modèles peuvent être un peu longs à fitter…).
moral1 <- brm(
response ~ 1,
data = d5,
family = cumulative("logit"),
cores = parallel::detectCores(),
)
moral2 <- brm(
response ~ 1 + action + intention + contact,
data = d5,
family = cumulative("logit"),
cores = parallel::detectCores()
)
Toutes les pentes sont négatives… ce qui signifie que chaque facteur réduit la réponse moyenne (i.e., le jugement de moralité). Ces pentes représentent des changements dans les log-odds cumulatifs.
moral1 <- add_criterion(moral1, "waic")
moral2 <- add_criterion(moral2, "waic")
w <- loo_compare(moral1, moral2, criterion = "waic")
w[,5:8] %>% round(digits = 2)
p_waic se_p_waic waic se_waic
moral2 9.13 0.04 37090.12 76.28
moral1 6.01 0.02 37854.46 57.67
model_weights(moral1, moral2, weights = "waic") %>%
round(digits = 2)
moral1 moral2
0 1
tidy(moral2, prob = 0.95)
term estimate std.error lower upper
1 b_Intercept[1] -2.839163e+00 0.04805513 -2.932892e+00 -2.744327e+00
2 b_Intercept[2] -2.156854e+00 0.04278314 -2.240424e+00 -2.074097e+00
3 b_Intercept[3] -1.574014e+00 0.03950732 -1.650146e+00 -1.498460e+00
4 b_Intercept[4] -5.524352e-01 0.03699212 -6.240745e-01 -4.799432e-01
5 b_Intercept[5] 1.164427e-01 0.03713198 4.442318e-02 1.895871e-01
6 b_Intercept[6] 1.023401e+00 0.03987168 9.436735e-01 1.100674e+00
7 b_action -7.093860e-01 0.04046395 -7.861238e-01 -6.308170e-01
8 b_intention -7.212895e-01 0.03688297 -7.940994e-01 -6.505540e-01
9 b_contact -9.614392e-01 0.05087764 -1.063520e+00 -8.625872e-01
10 lp__ -1.856178e+04 2.11921873 -1.856673e+04 -1.855863e+04
On peut représenter les prédictions du modèle en utilisant la fonction marginal_effects.
marg1 <- marginal_effects(moral2, "action", ordinal = TRUE)
p1 <-
plot(marg1, theme = theme_bw(base_size = 25), plot = FALSE)[[1]] +
scale_y_discrete(breaks = 1:7, labels = 1:7)
marg2 <- marginal_effects(moral2, "intention", ordinal = TRUE)
p2 <-
plot(marg2, theme = theme_bw(base_size = 25), plot = FALSE)[[1]] +
scale_y_discrete(breaks = 1:7, labels = 1:7)
marg3 <- marginal_effects(moral2, "contact", ordinal = TRUE)
p3 <-
plot(marg3, theme = theme_bw(base_size = 25), plot = FALSE)[[1]] +
scale_y_discrete(breaks = 1:7, labels = 1:7)
ggpubr::ggarrange(p1, p2, p3, nrow = 1, ncol = 3, common.legend = TRUE)
Posterior predictive checking. Pour plus d'informations sur la régression logistique ordinale, voir ce lien, ce tutorial, ou le chapitre 11 de McElreath (2015).
pp_check(moral2, nsamples = 1e2) + theme_bw(base_size = 25)